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A  LOCAL  DISCONTINUOUS  GALERKIN  METHOD  FOR  KDV-TYPE  EQUATIONS* 

JUE  YANt  AND  CHI-WANG  SHU^ 

Abstract.  In  this  paper  we  develop  a  local  discontinuous  Galerkin  method  for  solving  KdV  type  equa¬ 
tions  containing  third  derivative  terms  in  one  and  two  space  dimensions.  The  method  is  based  on  the 
framework  of  the  discontinuous  Galerkin  method  for  conservation  laws  and  the  local  discontinuous  Galerkin 
method  for  viscous  equations  containing  second  derivatives,  however  the  guiding  principle  for  inter-cell  fluxes 
and  nonlinear  stability  is  new.  We  prove  stability  and  a  cell  entropy  inequality  for  the  square  entropy 
for  a  class  of  nonlinear  PDEs  of  this  type  both  in  one  and  multiple  spatial  dimensions,  and  give  an  error 
estimate  for  the  linear  cases  in  the  one  dimensional  case.  The  stability  result  holds  in  the  limit  case  when  the 
coefficients  to  the  third  derivative  terms  vanish,  hence  the  method  is  especially  suitable  for  problems  which 
are  “convection  dominate” ,  i.e.  those  with  small  second  and  third  derivative  terms.  Numerical  examples  are 
shown  to  illustrate  the  capability  of  this  method.  The  method  has  the  usual  advantage  of  local  discontinuous 
Galerkin  methods,  namely  it  is  extremely  local  and  hence  efficient  for  parallel  implementations  and  easy  for 
h-p  adaptivity. 

Key  words,  discontinuous  Galerkin  method,  KdV  equation,  stability,  error  estimate 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  In  this  paper  we  develop  a  local  discontinuous  Galerkin  method  for  solving  KdV 
type  equations  containing  third  derivative  terms  in  one  and  multiple  spatial  dimensions.  An  example  of  such 
PDEs  is  the  original  KdV  equation  [18] 

ut  +  {au  + I3u^)x+(juxxx=^,  (1-1) 

where  a,  j3  and  a  are  constants.  Our  scheme  can  be  designed  and  proven  stable  for  more  general  nonlinear¬ 
ities,  namely 


ut  +  f{u)x  +  {r'{u)g{r{u)x)x)x  =  0  (1.2) 

in  one  space  dimension  for  arbitrary  (smooth)  functions  /,  g  and  r,  where  r'{u)  =  and 

d  d  /  d 

fi{u)xi  +  ri{u)'^gij{ri{u)^,)^. 

i^i  \  i=i 

in  multiple  spatial  dimensions  for  arbitrary  (smooth)  functions  /,,  gij  and  r,. 

KdV  type  equations  describe  the  propagation  of  waves  in  a  variety  of  nonlinear,  dispersive  media  and 
appear  often  in  applications.  See,  e.g.  [1].  Various  numerical  methods  have  been  proposed  and  used 
in  practice  to  solve  this  type  of  equations,  see,  e.g.  [4,  5,  17].  However,  in  many  situations,  such  as  in 
the  quantum  hydrodynamic  models  of  semiconductor  device  simulations  [15]  and  in  the  dispersive  limit  of 
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conservation  laws  [19],  the  third  derivative  terms  might  have  small  or  even  zero  coefficients  in  some  parts 
of  the  domain.  We  will  call  such  cases  as  “convection  dominated”.  The  design  of  stable,  efficient  and  high 
order  methods,  especially  those  for  the  “convection  dominated”  cases,  i.e.  when  the  third  derivative  terms 
are  small  (|(t|  <C  1  in  (1.1)),  remains  a  challenge. 

The  discontinuous  Galerkin  method  is  a  class  of  finite  element  methods  using  completely  discontinuous 
piecewise  polynomial  space  for  the  numerical  solution  and  the  test  functions.  One  certainly  needs  to  use 
more  degrees  of  freedom  because  of  the  discontinuities  at  the  element  boundaries,  however  this  also  gives  one 
a  room  to  design  suitable  inter-element  boundary  treatments  (the  so-called  fluxes)  to  obtain  highly  accurate 
and  stable  methods  in  many  difficult  situations. 

The  first  discontinuous  Galerkin  method  was  introduced  in  1973  by  Reed  and  Hill  [20],  in  the  framework 
of  neutron  transport  (steady  state  linear  hyperbolic  equations).  A  major  development  of  the  discontinuous 
Galerkin  method  was  carried  out  by  Cockburn  et  al.  in  a  series  of  papers  [10,  9,  7,  11],  in  which  they 
established  a  framework  to  easily  solve  nonlinear  time  dependent  hyperbolic  conservation  laws  (i.e.  (1.2) 
and  (1.3)  without  the  third  derivative  terms)  using  explicit,  nonlinear ly  stable  high  order  Runge-Kutta  time 
discretizations  [22]  and  discontinuous  Galerkin  discretization  in  space  with  exact  or  approximate  Riemann 
solvers  as  interface  fluxes  and  TVB  (total  variation  bounded)  nonlinear  limiters  [21]  to  achieve  non-oscillatory 
properties  for  strong  shocks. 

The  discontinuous  Galerkin  method  has  found  rapid  applications  in  such  diverse  areas  as  aeroacoustics, 
electro-magnetism,  gas  dynamics,  granular  flows,  magneto-hydrodynamics,  meteorology,  modeling  of  shallow 
water,  oceanography,  oil  recovery  simulation,  semiconductor  device  simulation,  transport  of  contaminant  in 
porous  media,  turbomachinery,  turbulent  flows,  viscoelastic  flows  and  weather  forecasting,  among  many 
others.  Good  references  for  the  discontinuous  Galerkin  method  and  its  recent  development  include  the 
survey  paper  [8],  other  papers  in  that  Springer  volume,  and  the  review  paper  [13]. 

The  original  discontinuous  Galerkin  method  was  designed  to  solve  first  order  hyperbolic  problems.  A 
simple  example  to  illustrate  the  essential  ideas  is  the  linear  transport  equation 

Ut+Ux=Q.  (1.4) 

Let’s  denote  the  mesh  by  Ij  =[xj_i,  for  j  =  with  the  center  of  the  cell  denoted  by  xj  = 

5  each  cell  by  Axj  =  Xj_^i  —  Xj_i.  We  will  denote  Ax  =  maxj  Axj.  If  we 

multiply  (1.4)  by  an  arbitrary  test  function  v(x),  integrate  over  the  interval  Ij,  and  integrate  by  parts,  we 
get 

/  Ufvdx  —  /  uvxdx  +  u{xj_^_i,t)v{xj_^_i)—u{xj_i,t)v{xj_i)=0.  (1-5) 

This  is  the  starting  point  for  designing  the  discontinuous  Galerkin  method.  We  replace  both  the  solution  u 
and  the  test  function  v  by  piecewise  polynomials  of  degree  at  most  k.  That  is,  u,v  £  Vax  where 

Vax  =  {u  :  u  is  a  polynomial  of  degree  at  most  A:  for  a;  €  Ij,  j  =  1, ...,  77}  .  (1.6) 

With  this  choice,  there  is  an  ambiguity  in  (1.5)  in  the  last  two  terms  involving  the  boundary  values  at  Xjj_^, 

as  both  the  solution  u  and  the  test  function  v  are  discontinuous  exactly  at  these  boundary  points.  The 
idea  is  to  treat  these  terms  by  an  upwinding  mechanism  (information  from  characteristics),  borrowed  from 
successful  high  resolution  finite  volume  schemes.  Thus  u  at  the  interfaces  Xjj_i  is  given  by  a  single  valued 
numerical  flux  Ujj_i  =  determined  from  upwinding,  and  v  at  the  interfaces  Xjj_^  are  given  by  the 
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values  taken  from  inside  the  cell  Ij,  namely  v  ^  and  Notice  that  we  use  v  and  to  denote  the  left 

and  right  limits  of  n,  respectively,  at  the  interface  where  v  is  discontinuous.  For  more  general  nonlinear  fluxes 
/(u),  the  only  difference  is  that  the  single  valued  flux  /j_,_  i  would  be  taken  as  a  monotone  flux  depending  on 
both  u~  i  and  on  ut,  i  (exact  or  approximate  Riemann  solvers  in  the  system  case).  The  resulting  method 

3\  2  3\  2 

of  the  lines  ODE  is  then  discretized  by  the  nonlinearly  stable  high  order  Runge-Kutta  time  discretizations 
[22].  Nonlinear  TVB  limiters  [21]  may  be  used  if  the  solution  contains  strong  discontinuities.  The  schemes 
thus  obtained,  for  solving  hyperbolic  conservation  laws  ((1.2)  and  (1.3)  without  the  third  derivative  terms), 
have  the  following  attractive  properties: 

1.  It  can  be  easily  designed  for  any  order  of  accuracy.  In  fact,  the  order  of  accuracy  can  be  locally 
determined  in  each  cell,  thus  allowing  for  efficient  p  adaptivity. 

2.  It  can  be  used  on  arbitrary  triangulations,  even  those  with  hanging  nodes,  thus  allowing  for  efficient 
h  adaptivity. 

3.  It  is  extremely  local  in  data  communications.  The  evolution  of  the  solution  in  each  cell  needs  to 
communicate  only  with  the  immediate  neighbors,  regardless  of  the  order  of  accuracy,  thus  allowing 
for  efficient  parallel  implementations.  See,  e.g.  [3]. 

4.  It  has  excellent  provable  nonlinear  stability.  One  can  prove  a  strong  stability  and  a  cell  entropy 
inequality  for  the  square  entropy,  for  the  general  nonlinear  cases,  for  any  orders  of  accuracy  on 
arbitrary  triangulations  in  any  space  dimension,  without  the  need  for  nonlinear  limiters  [16]. 

In  [12]  these  discontinuous  Galerkin  methods  were  generalized  to  solve  convection  diffusion  problems 
containing  second  derivative  terms.  This  was  motivated  by  the  successful  numerical  experiments  of  Bassi 
and  Rebay  [2]  for  the  compressible  Navier-Stokes  equations.  The  idea  can  be  illustrated  with  the  simple 
heat  equation 

'^XX  —  0  (^•*^) 

which  we  rewrite  into  a  first  order  system 


Ut-qx=Q,  q-Ux=Q, 


(1.8) 


we  can  then  formally  use  the  same  discontinuous  Galerkin  method  for  the  convection  equation  to  solve  (1.8), 
resulting  in  the  following  scheme:  find  u,q£  Vax  such  that,  for  all  test  functions  v,w  £  Vax, 


j  utvdx  +  j  qv^dx  -  qj_^_iv._^i  +  q. 


iv+  1  =  0 
^  2 


qwdx  + 


uwxdx  —  u,- 1  1  w  , ,  1  +  It,-  1  w'^  1  =  0. 
1+2  1-2 


(1.9) 


However,  there  is  no  longer  a  upwinding  mechanism  or  characteristics  to  guide  the  design  of  the  fluxes  Mj+i 
and  .  The  crucial  part  in  designing  a  stable  and  accurate  algorithm  (1.9)  is  a  correct  design  of  these 
fluxes.  In  [12],  criteria  are  given  for  these  fluxes  to  guarantee  stability,  convergence  and  a  sub-optimal  error 
estimate  of  order  k  for  piecewise  polynomials  of  degree  k.  The  (most  natural)  central  fluxes 

^i+h  =  ^  (vi  ’  <ij+h  =  ^  fci  +  <+ J 


would  satisfy  these  criteria  and  give  a  scheme  which  is  indeed  sub-optimal  in  the  order  of  accuracy  for  odd 
k  (i.e.  the  accuracy  is  order  k  rather  than  the  expected  order  A:  -I-  1  for  odd  k).  This  deficiency  however  is 
easily  removed  by  going  to  a  clever  choice  of  fluxes,  proposed  in  [12]: 


D+5 


j+i 


^i+h=^th- 


(1.11) 
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i.e.  we  alternatively  take  the  left  and  right  limits  for  the  fluxes  in  u  and  q  (we  could  of  course  also  take 
the  pair  and  q~_^_i  as  the  fluxes).  Notice  that  the  evaluation  of  (1.11)  is  simpler  than  that  of  the 

central  fluxes  in  (1.10),  and  this  easily  generalizes  to  multi  space  dimensions  on  arbitrary  triangulations. 
The  accuracy  now  becomes  the  optimal  order  A:  +  1  for  both  even  and  odd  k. 

The  schemes  thus  designed  for  the  heat  equation  (1.7),  or  in  fact  for  the  most  general  multi  dimensional 
nonlinear  convection  diffusion  equations  (nonlinear  both  in  the  first  derivative  convection  part  and  the  second 
derivation  diffusion  part) ,  retain  all  the  four  nice  properties  listed  above  for  the  method  used  on  convection 
equations.  Moreover,  the  appearance  of  the  auxiliary  variable  q  is  superficial:  when  a  local  basis  is  chosen 
in  cell  Ij  then  q  is  eliminated  and  the  actual  scheme  for  u  takes  a  form  similar  to  that  for  convection  alone. 
This  is  a  big  advantage  of  the  scheme  over  the  traditional  “mixed  methods” ,  and  it  is  the  reason  that  the 
scheme  is  termed  local  discontinuous  Galerkin  method  in  [12].  Even  though  the  auxiliary  variable  q  can  be 
locally  eliminated,  it  does  approximate  the  derivative  of  the  solution  u  to  the  same  order  of  accuracy,  thus 
matching  the  advantage  of  traditional  “mixed  methods”  on  this. 

The  purpose  of  this  paper  is  to  develop  a  similar  local  discontinuous  Galerkin  (LDG)  method  for  the 
KdV  like  equations  (1.1),  (1.2)  and  (1.3)  containing  third  derivative  terms.  Our  objective  is  to  design  the 
method  to  retain  again  all  the  four  nice  properties  listed  above  for  the  method  used  on  convection  and 
convection-diffusion  equations,  plus  the  feature  that  the  method  is  local,  namely  the  auxiliary  variables 
introduced  to  approximate  the  first  and  second  derivatives  of  the  solution  could  be  locally  eliminated. 

We  will  give  a  “preview”  of  the  method  on  the  simple  linear  equation 

Ut  Uxxx  —  0  (^•^^) 

which  we  again  rewrite  into  a  first  order  system 

Ut+Px=0,  p-qx=0,  q-ux=0.  (1.13) 


We  can  then  formally  use  the  same  discontinuous  Galerkin  method  for  the  convection  equation  to  solve 
(1.13),  resulting  in  the  following  scheme:  find  u,p,q  €  Vax  such  that,  for  all  test  functions  v,w,z  €  Va*, 


[  Utvdx-  [  pvxdx+p^+^v-^,-p^_^v+_,=0, 
j  pwdx  +  j  qwxdx  -  i  =  0, 

)  qzdx+)  uZxdx-u,,.z-^.+u, 

Jl.  J,.  J+2  ■'2^2 


(1.14) 


However,  the  fluxes  Pj+^,  be  designed  based  on  different  guiding  principles  than  the 

first  order  convection  or  second  order  diffusion  cases.  The  crucial  part  in  designing  a  stable  and  accurate 
algorithm  (1.14)  is  again  a  correct  design  of  these  fluxes.  It  turns  out  that  the  simple  choices 


i+r 


(1.15) 


would  guarantee  stability  and  convergence,  as  can  be  proven  later  in  this  paper  and  also  clearly  seen  in  Table 
1.1,  which  contains  numerical  and  L°°  errors  and  orders  of  accuracy  for  the  computed  solution  u  for  the 
method  (1.14)  with  the  fluxes  (1.15)  solving  the  equation  (1.12)  with  an  initial  condition  u(a;,  0)  =  sin(a;) 
over  the  interval  [0,27r]  and  periodic  boundary  conditions,  at  t  =  1,  using  uniform  meshes. 
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Table  1.1 


ut  +  Uxxx  =  0.  u{x,0)  =  sin(a;).  Periodic  boundary  conditions.  and  L°°  errors.  Uniform  meshes  with  N  cells.  LOG 
methods  with  A:  =  0, 1,  2,  3.  t  =  1. 


k 

N=10 

o 

II 

O 

II 

O 

00 

II 

error 

error 

order 

error 

order 

error 

order 

0 

2.2534E-01 

1.2042E-01 

0.91 

6.2185E-02 

0.95 

3.1582E-02 

0.98 

LOO 

4.3137E-01 

2.1977E-01 

0.97 

1.1082E-01 

0.98 

5.5376E-02 

1.00 

1 

1.7150E-02 

4.2865E-03 

2.00 

1.0716E-03 

2.00 

2.6792E-04 

1.99 

LOO 

5.8467E-02 

1.5757E-02 

1.89 

4.0487E-03 

1.96 

1.0210E-03 

1.99 

2 

8.5803E-04 

1.0823E-04 

2.98 

1.3559E-05 

2.99 

1.6958E-06 

3.00 

LOO 

4.0673E-03 

5.1029E-04 

2.99 

6.4490E-05 

2.98 

8.0722E-06 

3.00 

3 

3.3463E-05 

2.1035E-06 

3.99 

1.3166E-07 

3.99 

8.2365E-09 

3.99 

LOO 

1.8185E-04 

1.1157E-05 

3.97 

7.2362E-07 

3.99 

4.5593E-08 

3.99 

We  remark  that  the  choice  for  the  fluxes  (1.15)  is  not  unique.  In  fact,  the  crucial  part  is  to  take  p  and 
it  from  opposite  sides  and  to  take  q  from  the  right.  Thus 

Pi+ 1  =  P~+ 1 ,  Qj+^=  Qj+ 1 ,  «i+  i  =  1 , 

would  also  work. 

The  organization  of  the  paper  is  as  follows.  In  section  2  we  describe  the  method  for  the  one  dimensional 
case,  and  prove  its  nonlinear  stability  and  a  cell  entropy  inequality,  as  well  as  an  error  estimate  for  the 
linear  case.  In  section  3  multiple  spatial  dimensional  case  is  considered,  where  the  nonlinear  stability  is  given 
for  the  general  triangulations.  In  section  4  we  provide  several  numerical  examples  to  illustrate  the  capability 
of  the  method.  Concluding  remarks  and  remarks  about  future  work  are  given  in  section  5. 

2.  The  LDG  method  for  the  one  dimensional  case.  In  this  section,  we  present  and  analyze  the 
LDG  method  for  the  following  one  dimensional  nonlinear  problem: 

ut  +  f{u)x  +  {r'{u)g{r{u)x)x)x  =0,  0  <  a;  <  1  (2.1) 

with  an  initial  condition 

u{x,0)  =  u^{x),  0  <  a;  <  1  (2.2) 

and  periodic  boundary  conditions.  Here  f(u),  r(u)  and  gig)  are  arbitrary  (smooth)  nonlinear  functions. 
Notice  that  the  assumption  of  periodic  boundary  conditions  is  for  simplicity  only  and  is  not  essential:  the 
method  can  be  easily  designed  for  non-periodic  boundary  conditions.  Also  notice  that  the  linear  equation 
(1.12)  and  the  KdV  equation  (1.1)  are  both  special  cases  of  (2.1). 

To  define  the  LDG  method,  we  first  introduce  the  new  variables 

q  =  r{u)x,  P  =  g{q)x  (2.3) 

and  rewrite  the  equation  (2.1)  as  a  first  order  system: 

ut  +  {f{u)+r'{u)p)x=0,  p-g{q)x=0,  q-r{u)x=0.  (2.4) 
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The  LDG  method  is  obtained  by  discretizing  the  above  system  with  the  discontinuous  Galerkin  method.  This 
is  achieved  by  multiplying  the  three  equations  in  (2.4)  by  three  test  functions  v,w,z  respectively,  integrate 
over  the  interval  Ij,  and  integrate  by  parts.  We  also  need  to  pay  special  attention  to  the  boundary  terms 
resulting  from  the  procedure  of  integration  by  parts,  as  mentioned  in  the  previous  section.  Thus  we  seek 
piecewise  polynomial  solutions  u,p,q  €  Vax,  where  Vax  is  defined  in  (1.6),  such  that  for  all  test  functions 
v,w,z  e  Vax  we  have,  for  ^  <  j  <  N , 


Ufvdx  — 


{f{u)  +  r'{u)p)vxdx  +  (/  +  1  ^i+i  “  1 

J  pwdx  +  j  g{q)wxdx  - 


=  0, 
=  0, 
=  0. 


(2.5) 


Notice  that  we  still  use  letters  without  a  subscript  Aa;  to  denote  functions  in  the  finite  element  space  Vax,  to 
simplify  the  notations.  The  only  ambiguity  in  the  algorithm  (2.5)  now  is  the  definition  of  the  numerical  fluxes 
(the  “hats”),  which  should  be  designed  based  on  different  guiding  principles  than  the  first  order  convection 
or  second  order  diffusion  cases  to  ensure  stability.  It  turns  out  that  we  can  take  the  simple  choices  (we  omit 
the  subscripts  j  ±  ^  in  the  definition  of  the  fluxes  as  all  quantities  are  evaluated  at  the  interfaces 


f  =  f{u  ,u+), 


r(u~^)  —  r(u  )  ^  , 

-  P  =  P  , 


/.+  —  ?/. 


9  =  9{q  ),  r  =  r{u  ) 


(2.6) 


where  /(u“,u+)  is  a  monotone  flux  for  f{u),  namely  /(u“,u+)  is  a  Lipschitz  continuous  function  in  both 
arguments  u~  and  is  consistent  with  f{u)  in  the  sense  that  f{u,u)  =  f{u),  and  is  a  non-decreasing 
function  in  u~  and  a  non-increasing  function  in  u~^ .  Likewise,  —g(q~,q~^)  is  a  monotone  flux  for  —g(q), 
namely  g{q~,q~^)  is  a  Lipschitz  continuous  function  in  both  arguments  q~  and  is  consistent  with  g{q) 
in  the  sense  that  g{q,  q)  =  g{q),  and  is  a  non-increasing  function  in  q~  and  a  non-decreasing  function  in  g+. 
Examples  of  monotone  fluxes  which  are  suitable  for  discontinuous  Galerkin  methods  can  be  found  in,  e.g. 
[10].  We  could  for  example  use  the  simple  Lax-Friedrichs  flux 


f{u  ,u+)  =  -(/(u  )  +  f{u'+)  -  a{u~^ -u  ))  ,  a  =  max|/'(u)|. 


(2.7) 


where  the  maximum  is  taken  over  a  relevant  range  of  u.  The  algorithm  is  now  well  defined. 

We  remark  that  the  choice  for  the  fluxes  (2.6)  is  not  unique.  In  fact,  the  crucial  part  is  to  take  p  and  f 
from  opposite  sides.  Thus 


f  =  f{u  ,u+), 


-  r{u+)-r{u  )  .  _  . 

r' = - - - ,  p  =  p  ,  g  =  g{q  ,q^),  r  =  r{u^) 

?/.“•"  —  ?/. 


would  also  work. 

We  also  remark  that  the  algorithm  (2.5)-(2.6)  is  very  easy  for  numerical  implementation.  Given  u,  one 
first  uses  the  third  equation  in  (2.5)  to  obtain  q.  This  is  achieved  locally:  q  in  Ij  can  be  obtained  with  the 
information  of  u  in  the  cells  Ij  and  Ij-i-  The  second  equation  in  (2.5)  is  then  used  to  obtain  p  locally:  p  in 
Ij  can  be  obtained  with  the  information  of  q  in  (at  most)  the  cells  Ij,  Ij-i  and  /j+i.  Finally,  the  update  of 
the  solution  u  is  obtained  using  the  first  equation  in  (2.5),  again  locally,  namely  the  update  of  u  in  Ij  can 
be  obtained  with  the  information  of  u  in  (at  most)  the  cells  Ij,  Ij-i  and  /j+i  and  that  of  p  in  the  cells  Ij 
and  Ij+i . 
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We  have  the  following  “cell  entropy  inequality”  for  the  scheme  (2.5)-(2.6).  This  is  a  generalization  of 
the  cell  entropy  inequality  obtained  in  [16]  for  the  discontinuous  Galerkin  method  applied  to  hyperbolic 
conservation  laws  (equation  (2.1)  with  g{q)  =  0). 


Proposition  2.1.  (cell  entropy  inequality)  There  exist  numerical  entropy  fluxes  such  that  the  solution 

to  the  scheme  (2.5)-(2.6)  satisfies 


£ 

dt 


) 


2  )  — 


<  0. 


Proof:  We  sum  up  the  three  equalities  in  (2.5)  and  introduce  the  notation 

=  f  -  /(/(„)  +  /(„),).,*  +  (/  + 


pwdx  +  L  g{q)w^dx  -  g^j^i 


■2  i+| 

+  I  d^dx  +  r(u)z,dx  -  ■ 

Clearly,  the  solutions  u,  p,  q  of  the  scheme  (2.5)-(2.6)  satisfy 

Bj{u,p,q-,v,w,z)  =  0 

for  all  v,w,z  €  Vax-  We  then  take 


(2.8) 


(2.9) 


(2.10) 


V  =  w  =  z  =  —p 

to  obtain,  after  some  algebraic  manipulations, 

d  f  f  ii^ix  \  /  •'  \ 

0  =  Bj{u,p,q-,u,q,  -p)  = -^  J ^  dx  +  [Hj+i  -  Hj_ij  +  0j_i 

with  the  numerical  entropy  flux  H  defined  by 

H  =  —F{u~)  +  G{q~)  —  r{u~)p~  +  ^/  +  r'pj  u~  —  gq~  +  rp~ 
and  the  extra  term  0  given  by 

0  =  [F(u)  -  G(q)  +  r(u)p]  -  (^f  +  r'p^  [u]  +  g[q]  -  f[p], 

Here 

Fiu)  =  j  f{u)du,  G{q)  =  j  g{q)dq, 

and 


[u]  =  v~^  —  V 

denotes  the  jump  of  v.  Notice  that  we  have  dropped  the  subscripts  about  the  location  j  —  ^  or  j  +  i  as 
all  these  quantities  are  defined  at  a  single  interface  and  depend  only  on  the  left  and  right  values  at  that 
interface.  Now  all  we  need  to  do  is  to  verify  0  >  0.  To  this  end,  we  notice  that,  with  the  definition  (2.6)  of 
the  numerical  fluxes  and  with  simple  algebraic  manipulations,  we  easily  obtain 

[r(u)p]  —  r'p[u]  —  f[p]  =  0 
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and  hence 


0  =  [F(u)]  -  f[u]  -  [G(g)]  + : 


[  (f(s)-f(u  ,u+)]ds-  f  {g(s)-g(q  ,g+)) 

J U~  ^  '  J Q~ 


ds 


(2.11) 


>0, 


where  the  last  inequality  follows  from  the  monotonicity  of  the  fluxes  /  and  —g.  This  finishes  the  proof.  □ 

Now  the  stability  of  the  method  is  a  trivial  corollary: 

Corollary  2.1.  (L^  stability)  The  solution  to  the  scheme  (2.5)-(2.6)  satisfies  the  stability 

d  /  u‘^{x,t) 


dt  Jo 


dx  <  0. 


(2.12) 


Proof:  We  simply  add  up  (2.8)  over  j.  □ 

About  time  discretizations,  if  we  denote  the  semi-discrete  LDG  method  (2.5)-(2.6)  by 

ut  =  R(u), 


then  the  following  implicit  0  scheme 

„.n+l  _  „.n 

- — - =i?(u”+®),  =  (1  -  6»)u”  +  6»u”+i  (2.13) 

will  also  satisfy  the  same  cell  entropy  inequality  and  stability  as  long  as  ^  <  0  <  1-  Notice  that  this 
includes  the  first  order  backward  Euler  and  second  order  Crank-Nicolson  implicit  time  discretizations  as 
special  cases.  See  [16]  for  the  purely  hyperbolic  case. 


Proposition  2.2.  (implicit  time  discretization)  The  cell  entropy  inequality  and  the  stability  also  hold 
for  the  time  discretization  (2.13)  with  i  <  0  <  1  for  the  scheme  (2.5)-(2.6).  That  is, 


(u"+^(a;))^  —  (u"(x)y 
2At 


dx  + 

J+5 


_ 

J~2 


<  0, 


(2.14) 


and 


f\u"+^(x)fdx  <  f\u"{x)ydx. 

Jo  Jo 


(2.15) 


Proof:  If  we  take  the  test  functions  at  n  +  e.g.  v  =  given  by  (2.13),  we  obtain  just  as  before 

f  i,n+.  _  ^.+,  < 

/  r.  At  2  JJ  2 


which  can  be  rewritten  as 

'  (u"+^(a;))^  —  (u"(x))‘^ 
2At 


dx  +  m+l-m+l  +  [0-- 


1 


(u"+^(a;)  —  u"(a;))^ 
At 


dx  <  0. 


Thus,  a  sufficient  condition  to  get  the  cell  entropy  inequality  (2.14)  is  just  0  >  ^-  Again,  (2.15)  is 
obtained  simply  by  adding  up  (2.14)  over  j.  □ 


The  stability  result  obtained  here  can  be  used  to  get  an  error  estimate  in  for  the  numerical  solution  u, 
when  the  equation  (2.1)  is  linear.  Without  loss  of  generality  we  may  take  f{u)  =  u,  g{q)  =  q  and  r{u)  =  u, 
resulting  in  the  equation 

<+<+<..  =  0.  (2.16) 

Noticed  that  we  have  used  the  notation  to  denote  the  exact  solution  of  the  PDE  in  order  not  to  confuse 
with  the  numerical  solutions.  We  have  the  following  result,  where  C  here  and  below  denotes  a  generic 
constant  which  may  be  of  different  values  at  different  locations. 

Proposition  2.3.  (error  estimate)  The  error  for  the  scheme  (2.5)-(2.6)  applied  to  the  linear  PDE  (2.16) 
satisfies 

<C'Aa;*+^  (2.17) 

where  the  constant  C  depends  on  the  derivatives  of  and  time  t. 

Proof:  Eirst,  we  notice  that,  in  this  linear  case,  most  monotone  fluxes  simply  become  upwinding 

/(u“,u+)=u“,  g{q~  ,q'^)  =  q'^ , 

and  this  is  what  we  will  assume.  It  is  then  easy  to  work  out  the  exact  form  of  0  in  (2.11)  for  the  cell  entropy 
inequality  to  be 

0  =  +  [#)•  (2.18) 

We  now  notice  that  the  exact  solution  of  the  PDE  (2.16),  u'^ ,  q'^  =  u%  and  p'^  =  clearly  satisfies 

Bj{u^  ,p‘^  ,q^-,v,w,z)  =  0 

for  all  v,w,z  €  Vax,  where  Bj  is  defined  by  (2.9).  Taking  the  difference  between  the  above  equality  and 
(2.10),  we  obtain  the  error  equation 

Bj{u^  —  u,p^  —  p,  q'^  —  q;  v,w,  z)  =  0  (2-19) 

for  all  v,w,z  €  Vax-  As  usual  this  error  equation  is  the  basic  starting  point  of  error  estimates. 

We  now  take 


v  =  Su^-u,  w  =  Vq^-q,  z=p-Vp\  (2.20) 

in  the  error  equation  (2.19).  Here  V  is  the  standard  projection  into  Va*,  that  is,  for  each  j, 

/  (Vw(x)  —  w(x))v(x)dx  =  0  Vn  e  P*, 

Jii 

where  P*  denotes  the  space  of  all  polynomials  of  degree  at  most  k.  In  other  words,  the  difference  between 
the  projection  Vw  and  the  original  function  w  is  orthogonal  to  all  polynomials  of  degree  up  to  k  in  each 
interval.  5  is  a  special  projection  into  Va*  which  satisfies,  for  each  j, 


j  {Sw{x)  —  w{x))v{x)dx  =  \lv  &  ^  and 
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in  other  words,  the  difference  between  the  projection  Sw  and  the  original  function  w  is  orthogonal  to  all 
polynomials  of  degree  up  to  A:  —  1  in  each  interval,  and  the  projection  agrees  with  the  function  at  the  right 
boundary  in  each  interval.  This  special  projection  is  needed  for  u  because  we  have  no  control  on  the  jumps 
of  p  in  the  cell  entropy  inequality,  see  (2.18).  Substituting  (2.20)  into  the  error  equation  (2.19)  and  moving 
terms,  we  obtain 


Bj(v,  —z,w;v,w,  z)  =  Bj(v'^ ,  —z'^,w'^;  v,w,z) 


(2.21) 


where  v,  w,  z  are  given  by  (2.20),  and  ,  z^  are  given  by 

v^=Su^-u\  z^=p^-Vp^.  (2.22) 

By  the  same  argument  as  that  used  for  the  cell  entropy  inequality,  the  left  hand  side  of  (2.21)  becomes 

Bj{v,-z,w,v,w,z)  =  -^  I  dx+  (^1 

where,  by  (2.18), 

%-i  =  I  (m3_j  +  [»i3_j) . 

The  right  hand  side  of  (2.21)  can  be  written  out  as 

Bj(u",  -.z", w" ■,v,w,z)  =I  +  II  +  TTI  +  TV 

where 


(2.23) 


(2.24) 


(2.25) 


T  =  /  VfVdx, 


(2.26) 


XT  =  —  /  z^wdx+  /  zdx  —  /  (v'^  —  z'^)vxdx  +  /  w'^Wxdx  +  /  v^Zxdx, 
Jij  Jii  Jii  Jij  Jii 

i^u)  -  K-i + YuY  + Y-i) 


III = - 


(2.27) 


(2.28) 


and 


Tv  =  y+i-y_i  (2.29) 

for  some  flux  function  h.  Notice  that  v,w,z  are  given  by  (2.20)  and  ,x^  are  given  by  (2.22),  respectively. 

Now,  by  using  the  simple  inequality  ah  <  \{a^  +  6^),  and  standard  approximation  theory  on  = 
{Su'^  —  see,  e.g.  [6],  we  have 

dx. 

Because  P  is  a  local  projection,  and  S,  even  though  not  a  local  projection,  does  have  the  property  that 
w  —  Sw  is  locally  orthogonal  to  all  polynomials  of  degree  up  to  A:  — 1,  all  the  terms  mXX  are  actually  zero.  The 
last  term  in  XXX  is  zero,  because  of  the  special  interpolating  property  of  the  projection  S.  An  application 


X<CAxf+^^  +  1^ 
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of  the  simple  inequality  ab  <  +b‘^)  for  the  first  two  terms  in  XXX  and  standard  approximation  theory 

on  the  point  values  of  =  {Su^  —  u^)  +  {Vp^  —  p^)  and  of  =  Vq^  —  (see,  e.g.  [6])  then  gives 

XXX  <  CiAxfJy^  +  Axf+‘^  +  ^  ([u]2  +  H^)  . 

Finally,  XV  only  contains  flux  difference  terms  which  will  vanish  upon  a  summation  in  j. 

Combining  all  these  and  summing  over  j  we  obtain  the  following  inequality 

An  integration  in  t  plus  the  standard  approximation  theory  on  v'^  =  Su'^  —  u'^  then  gives  the  desired  error 
estimate  (2.17).  □ 

3.  The  LDG  method  for  the  multiple  dimensional  case.  In  this  section,  we  generalize  the  scheme 
discussed  in  the  previous  section  to  multiple  spatial  dimensions  x  =  {xi,  -  ■  ■  , Xd) .  We  solve  the  following 
nonlinear  problem: 

d  d  /  d 

fi{u)xi  +  Y1  r-{u)'^gij{ri{u)^,)^. 
i^i  y  i=i 

with  an  initial  condition 


Ut  +  Y. 


=  0,  0  <  a;,  <  1, 


=  1, 


(3.1) 


u{x,  0)  =  u^{x),  0  <  a;,  <  1,  i  =  1,  •  •  •  ,d  (3-2) 

and  periodic  boundary  conditions.  Here  fi(u),  ri(u)  and  gij(q)  are  arbitrary  (smooth)  nonlinear  functions. 
Notice  that  the  assumption  of  a  box  geometry  and  periodic  boundary  conditions  is  for  simplicity  only  and 
is  not  essential:  the  method  can  be  easily  designed  for  arbitrary  domain  and  for  non-periodic  boundary 
conditions. 

Let’s  denote  a  triangulation  of  the  unit  box  by  Tax,  consisting  of  non-overlapping  polyhedra  covering 
completely  the  unit  box.  Hanging  nodes  are  allowed.  Here  Ax  measures  the  longest  edge  of  all  polyhedra  in 
Tax-  We  again  denote  the  finite  element  space  by 

^Ax  =  {v  ■  u  is  a  polynomial  of  degree  at  most  A:  for  a;  €  K,  \fK  €  Tax}  ■  (3.3) 

Similar  to  the  one  dimensional  case,  to  define  the  LDG  method,  we  first  introduce  the  new  variables 

d 

qi=ri{u)xi,  P*  =  i  =  (3.4) 

i=i 

and  rewrite  the  equation  (3.1)  as  a  first  order  system: 

d 

ut  +  +  ‘r'M)Pi)xi  =  0, 

d 

P*  -  =  0,  qi  -  ri{u)xi  =  0,  i  =  l,---,d.  (3.5) 

i=i 

The  LDG  method  is  obtained  by  discretizing  the  above  system  with  the  discontinuous  Galerkin  method. 
This  is  achieved  by  multiplying  the  equations  in  (3.5)  by  test  functions  v,Wi,Zi  respectively,  integrate  over 
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an  element  K  €  Ta*,  and  integrate  by  parts.  We  again  need  to  pay  special  attention  to  the  boundary 
terms  resulting  from  the  procedure  of  integration  by  parts,  as  in  the  one  dimensional  case.  Thus  we  seek 
piecewise  polynomial  solutions  u,pi,qi  €  where  is  defined  in  (3.3),  such  that  for  all  test  functions 
v,Wi,Zi  e  we  have 


f  utvdx-'^f  {fi{u)+r[{u)pi)vxidx+f  ds  =  0, 

Jk  Jk  JdK 

/  piWidx  +  '^  /  gijiqi)iwi)x^dx  -  /  ^ 

Jk  Jk  JdK 

/  ri(u)(zi) 

Jk 


=  0, 


L 


i=i 
qiZidx  + 


,dx  — 


zi^tKds  =  0, 


=  I,  -  ,d 


1 ,  •  •  •  ,  d. 


(3.6) 


where  dK  is  the  boundary  of  element  K,  and  the  numerical  fluxes  (the  “hats”)  are  defined  similar  to  the 
one  dimensional  cases,  namely 

rii^x 


h„^  =/„^. if  + 


gi,nK  =  9i,nK,K{q  ,q^  ’"),  ri^nK=ri(u  )ni^K- 


(3.7) 


Here  tik  =  (iii.if,  •  •  •  ,nd^K)  is  the  outward  unit  normal  for  element  K  along  the  element  boundary  dK, 
^tntK  denotes  the  value  of  u  evaluated  from  inside  the  element  K,  and  denotes  the  value  of  u  evaluated 
from  outside  the  element  K  (inside  the  neighboring  element).  On  the  other  hand,  denotes  the  value  of 
p  evaluated  from  a  pre-designated  “plus”  side  along  an  edge  e,  which  is  always  the  boundaries  of  two 
neighboring  elements.  For  example,  we  could  choose  a  fixed  vector  jj  which  is  not  parallel  with  any  normals 
of  element  boundaries,  and  then  designate  the  “plus”  side  to  be  the  side  at  the  end  of  the  arrow  of  the  normal 
n  with  n-/?  >  0,  see  Figure  3.1.  fuK ,k(u^"^‘^  is  a  monotone  flux  for  /„^  (u)  =  fi(u)ni^K,  namely 

fnK,K(u^"^‘^ ,  is  a  Lipschitz  continuous  function  in  both  arguments  and  ,  is  consistent  with 

/„^  (u)  in  the  sense  that  /„^  {u,u)  =  /„^  (u),  and  is  a  non-decreasing  function  in  and  a  non-increasing 

function  in  u*^**^.  Moreover,  it  is  conservative  (that  is,  there  is  only  one  flux  at  each  edge  shared  by  two 
elements,  added  to  the  residue  for  one  and  subtracted  from  the  reside  for  another),  namely 

fnK,K{a,b)  =  -/„^,,if'(6,a) 

where  K  and  K'  share  the  same  edge  where  the  flux  is  computed  and  hence  hk'  =  —uk-  Likewise, 
-gXn^K<yqT'*’^  ,qT*’^)  is  a  monotone  flux  for  -5f*,„^(g*)  =  -J2j=i9ij(g)nj,K-  Notice  that  we  can  again 
use  the  one  dimensional  monotone  fluxes  as  in  the  previous  section.  For  example,  we  can  use  the  simple 
Lax-Friedrichs  flux 

^  ,  (3.8) 

a  =  max|/;  (u)|, 

U 

where  the  maximum  is  taken  over  a  relevant  range  of  u.  The  algorithm  is  now  well  defined. 

Again,  the  algorithm  (3.6)-(3.7)  is  very  easy  for  numerical  implementation.  Given  u,  one  first  locally 
solves  for  the  g,,  then  locally  solves  for  the  pi,  and  finally  locally  solves  for  the  update  of  u.  All  the  advantages 
listed  for  the  method  for  the  one  dimensional  case  are  still  valid  in  this  multiple  dimensional  case. 

We  still  have  the  following  “cell  entropy  inequality”  for  the  scheme  (3.6)-(3.7).  The  proof  follows  the 
same  lines  as  that  for  the  one  dimensional  case,  so  we  will  omit  it. 


12 


p 


T 


Fig.  3.1.  Illustration  of  the  definition  of  “plus”  and  “minus”  sides  determined  by  a  pre- determined  veetor  /3. 

Proposition  3.1.  (cell  entropy  inequality)  There  exist  conservative  numerical  entropy  fluxes  such 

that  the  solution  to  the  scheme  (3.6)-(3.7)  satisfles 

□ 


The  stability  of  the  method  is  then  again  a  trivial  corollary,  by  summing  up  the  cell  entropy  inequal¬ 
ities  over  K: 


Corollary  3.1.  stability)  The  solution  to  the  scheme  (3.6)-(3.7)  satisfles  the  stability 


£ 

dt 


I 


\x,t) 


dx  <  0. 


(3.10) 


□ 


The  same  cell  entropy  inequality  also  holds  for  the  implicit  time  discretizations: 


Proposition  3.2.  (implicit  time  discretization)  The  cell  entropy  inequality  and  the  stability  also  hold 
for  the  time  discretization  (2.13)  with  i  <  0  <  1  for  the  scheme  (3.6)-(3.7).  That  is, 

'(u"+^(a;))^  —  (u"(x))‘^' 


K 


2At 


f 

JdK 


dx+l  Kt'Kds<Q, 


and 


□ 


[  (u^+Hx)fdx  <  [  (u^(x))^dx. 

Jn  Jn 


(3.11) 


(3.12) 


Unfortunately,  we  could  not  get  the  optimal  error  estimate  because  of  the  lack  of  a  suitable  projection 
S  similar  to  the  one  dimensional  case.  However,  numerical  examples  in  the  next  section  verify  that  the 
accuracy  holds  as  in  the  one  dimensional  case. 
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Table  4.1 


ut  +  Uxxx  =  0.  u{x,0)  =  sin(a;).  Periodic  boundary  conditions.  and  L°°  errors.  Non-uniform  meshes  with  N  cells. 
LDG  methods  with  A:  =  0, 1,  2,  3.  t  =  1. 


k 

N=10 

O 

II 

O 

II 

O 

00 

II 

error 

error 

order 

error 

order 

error 

order 

0 

2.2222E-01 

1.2014E-01 

0.88 

6.2532E-02 

0.94 

3.1900E-02 

0.97 

LOO 

4.3282E-01 

2.2006E-01 

0.97 

1.1210E-01 

0.97 

5.8810E-02 

0.93 

1 

2.0144E-02 

5.2347E-03 

1.94 

1.3322E-03 

1.97 

3.3592E-04 

1.98 

LOO 

8.8110E-02 

2.3302E-02 

1.93 

5.9387E-03 

1.97 

1.4969E-03 

1.98 

2 

9.8394E-04 

1.1974E-04 

3.03 

1.4953E-05 

3.00 

1.8687E-06 

3.00 

LOO 

5.2984E-03 

6.8421E-04 

2.95 

8.5138E-05 

3.00 

1.0728E-05 

2.99 

3 

7.3589E-05 

4.6509E-06 

3.98 

2.9191E-06 

3.99 

2.0141E-08 

3.86 

LOO 

3.4438E-04 

2.2260E-05 

3.95 

1.3992E-06 

3.99 

9.1039E-08 

3.94 

4.  Numerical  examples.  In  this  section  we  provide  a  few  preliminary  numerical  examples  to  illustrate 
the  accuracy  and  capability  of  the  method.  Attention  has  not  been  paid  to  efficiency  in  time  discretizations, 
so  explicit  third  order  Runge-Kutta  method  [22]  is  used.  Study  of  suitable  implicit  time  discretizations  which 
have  efficient  iterative  solvers  maintaining  the  local  structure  of  the  method  is  the  subject  of  a  future  study. 

We  would  like  to  illustrate  through  these  numerical  examples  the  high  order  accuracy  of  the  methods 
for  both  one  dimensional  and  two  dimensional,  both  linear  and  nonlinear  problems.  We  would  also  like  to 
illustrate  the  good  behavior  of  the  method  for  the  so-called  convection  dominated  cases,  namely  the  case 
where  the  coefficients  of  the  third  derivative  terms  are  small. 

Example  4.1.  We  compute  the  solution  of  the  linear  one  dimensional  equation 

Ut  Uxxx  —  0  (^•^) 

with  an  initial  condition  u{x,  0)  =  sin(a;)  and  periodic  boundary  conditions  (with  27r  periodicity).  The  exact 
solution  is  given  by  u{x,t)  =  sin(a;  +  t).  Both  uniform  meshes  and  non-uniform  meshes  are  used.  The 
non-uniform  meshes  in  this  and  later  examples  are  a  repeated  pattern  of  0.9Aa;  and  l.lAa;  with  an  even 
number  of  elements.  The  and  L°°  errors  and  the  numerical  order  of  accuracy  are  contained  in  Table  1.1 
(in  section  1)  for  the  uniform  mesh  case,  and  in  Table  4.1  for  the  non-uniform  mesh  case.  We  can  clearly 
see  that  the  method  with  P*  elements  are  giving  a  uniform  {k  l)-th  order  of  accuracy  in  both  norms  for 
both  the  uniform  and  the  non-uniform  meshes. 


Example  4.2.  We  compute  the  solution  of  the  linear  two  dimensional  equation 

Ut  T  Uxxx  T  '^yyy  —  0  (^•^) 

with  an  initial  condition  u{x,y,0)  =  sin(a;  -I-  y)  and  periodic  boundary  conditions  (with  27r  periodicity)  in 
both  directions.  The  exact  solution  is  given  by  u{x,y,t)  =  sin(a;  y  Both  uniform  and  non-uniform 

rectangular  meshes  are  used.  The  non-uniform  meshes  are  a  repeated  pattern  of  0.9Aa;  and  l.lAa;,  in  both 
directions,  with  an  even  number  of  edges  in  both  directions.  The  and  L°°  errors  and  the  numerical  order 
of  accuracy  are  contained  in  Table  4.2  for  the  uniform  mesh  case,  and  in  Table  4.3  for  the  non-uniform  mesh 
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Table  4.2 


ut  +  Uxxx  +  Uyyy  =  0.  u{x,y,0)  =  sin(a;  +  y).  Periodic  boundary  conditions.  and  L°°  errors.  Uniform  meshes  with 
N  X  N  cells.  LOG  methods  with  A:  =  0, 1,  2,  3.  t  =  1. 


k 

10x10 

20x20 

40x40 

error 

error 

order 

error 

order 

0 

3.5528E-01 

2.0535E-01 

0.79 

1.1090E-01 

0.89 

LOO 

7.1359E-01 

4.0190E-01 

0.82 

2.1165E-01 

0.92 

1 

3.3603E-02 

9.0904E-03 

1.89 

2.4084E-03 

1.92 

LOO 

2.2074E-01 

6.1899E-02 

1.83 

1.5962E-02 

1.95 

2 

3.8750E-03 

4.8463E-04 

2.99 

6.0501E-05 

3.00 

LOO 

3.9084E-02 

4.8902E-03 

2.99 

6.1274E-04 

2.99 

3 

4.1491E-04 

2.6426E-05 

3.97 

1.6550E-06 

3.99 

LOO 

4.2847E-03 

2.8478E-04 

3.91 

1.7846E-05 

3.99 

Table  4.3 

ut  +  Uxxx  +  Uyyy  =  0.  u{x,y,0)  =  sin(a;  +  2/).  Periodic  boundary  conditions.  and  L°°  errors.  Non-uniform  meshes 
with  N  X  N  cells.  LOG  methods  with  A:  =  0,1,  2, 3.  t  =  1. 


k 

10x10 

20x20 

40x40 

error 

error 

order 

error 

order 

0 

3.5963E-01 

2.0788E-01 

0.79 

1.1228E-01 

0.88 

LOO 

7.3869E-01 

4.0713E-01 

0.85 

2.1681E-01 

0.91 

1 

3.4590E-02 

9.1681E-03 

1.92 

2.3412E-03 

1.97 

LOO 

2.5815E-01 

7.2978E-02 

1.82 

1.8533E-02 

1.97 

2 

4.0949E-03 

5.1285E-04 

2.99 

6.4054E-05 

3.00 

LOO 

5.0429E-02 

6.3078E-03 

2.99 

8.0584E-04 

2.97 

3 

4.5434E-04 

2.8854E-05 

3.97 

1.8080E-06 

3.99 

LOO 

6.0982E-03 

4.0321E-04 

3.92 

2.5340E-05 

3.99 

case.  We  can  clearly  see  again  that  the  method  with  P*  elements  are  giving  a  uniform  {k  +  l)-th  order  of 
accuracy  for  both  the  uniform  and  the  non-uniform  meshes. 


Example  4.3.  In  order  to  see  the  accuracy  of  the  method  for  nonlinear  problems,  we  compute  the  classical 
soliton  solution  of  the  KdV  equation 


tit  ^  )a;  Uxxx  —  0  (4-3) 

in  — 10  <  a;  <  12.  The  initial  condition  is  given  by 

u(a;,  0)  =  — 2sech^(a;), 


The  exact  solution  is 


u{x,t)  =  —  2  sech^(a;  —  4t). 
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Table  4.4 

The  KdV  equation  ut  —  +  Uxxx  =  0.  u(x,0)  =  —2seeh^(x).  Boundary  eondition  (4-4)-  errors. 

Uniform  meshes  with  N  eells.  LOG  methods  with  A:  =  0,1,2,  3.  t=  0.5. 


k 

o 

II 

O 

00 

II 

N=160 

N=320 

error 

error 

order 

error 

order 

error 

order 

0 

2.5292E-01 

1.9098E-01 

0.40 

1.3019E-01 

0.55 

7.9780E-02 

0.71 

LOO 

9.0170E-01 

6.8651E-01 

0.39 

4.6405E-01 

0.56 

2.8531E-01 

0.70 

1 

2.6512E-02 

4.6652E-03 

2.50 

1.0108E-03 

2.20 

2.5906E-04 

1.96 

LOO 

1.4748E-01 

3.4625E-02 

2.09 

1.1840E-02 

1.55 

3.3239E-03 

1.83 

2 

1.5317E-03 

1.8083E-04 

3.08 

2.2642E-05 

2.99 

2.8335E-06 

2.99 

LOO 

1.7486E-02 

2.7505E-03 

2.66 

3.5575E-04 

2.95 

4.4397E-05 

3.00 

3 

2.0631E-04 

1.3981E-05 

3.88 

8.9054E-07 

3.97 

5.6029E-08 

3.99 

LOO 

2.0155E-03 

2.1462E-04 

3.23 

1.4461E-05 

3.89 

9.1140E-07 

3.98 

Table  4.5 

The  KdV  equation  ut  —  3(t(^)^  +  Uxxx  =  0.  u(x,0)  =  —2seeh^(x).  Boundary  eondition  (4-4)-  o,nd  L°°  errors. 
Non-uniform  meshes  with  N  eells.  LOG  methods  with  A:  =  0,1, 2,  3.  t  =  0.5. 


k 

o 

II 

O 

00 

II 

N=160 

N=320 

error 

error 

order 

error 

order 

error 

order 

0 

2.4530E-01 

1.9004E-01 

0.37 

1.3390E-01 

0.50 

8.4635E-02 

0.66 

LOO 

1.0172E-I-00 

7.6826E-01 

0.40 

5.3383E-01 

0.52 

3.3655E-01 

0.66 

1 

2.7042E-02 

4.9065E-03 

2.46 

1.0555E-03 

2.21 

2.6978E-04 

1.97 

LOO 

1.4490E-01 

4.1570E-02 

1.80 

1.3925E-02 

1.57 

3.9129E-03 

1.83 

2 

1.9493E-03 

2.0134E-04 

3.27 

2.4926E-05 

3.01 

3.1208E-06 

2.99 

LOO 

2.2876E-02 

3.5163E-03 

2.70 

4.7161E-04 

2.89 

5.9033E-05 

2.99 

3 

3.0402E-04 

1.5462E-05 

4.29 

1.0064E-06 

3.94 

6.3370E-08 

3.99 

LOO 

2.7735E-03 

2.1464E-04 

3.69 

1.8358E-05 

3.55 

1.3119E-06 

3.80 

We  compute  the  solution  with  two  different  boundary  conditions.  Table  4.4  (uniform  mesh)  and  Table  4.5 
(non-uniform  mesh)  give  the  errors  of  numerical  solution  at  t  =  0.5  using  the  boundary  condition 

u{-10,t)  =  gi{t),  Ux{12,t)  =  g2{t),  Uxx{'i-2,t)  =  gsit)  (4.4) 

where  gi{t)  corresponds  to  the  data  from  the  exact  solution.  Notice  that  the  LDG  method  allows  an  easy 
implementation  of  such  boundary  conditions.  Table  4.6  (uniform  mesh)  and  Table  4.7  (non-uniform  mesh) 
give  the  errors  of  numerical  solution  using  the  periodic  boundary  conditions.  Although  the  exact  solution  is 
not  periodic,  the  large  size  of  the  computational  domain  allows  the  usage  of  periodic  boundary  conditions 
with  negligible  error.  We  can  see  from  these  tables  that  the  orders  of  accuracy  are  comparable  to  that  for 
the  linear  case. 
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Table  4.6 

The  KdV  equation  ut  —  3  (u^)^  +  Uxxx  =  0.  u(x,0)  =  —2seeh^(x).  Periodie  boundary  eondition.  and  L°°  errors. 
Uniform  meshes  with  N  eells.  LOG  methods  with  A:  =  0,1,2,  3.  t=  0.5. 


k 

o 

II 

O 

00 

II 

N=160 

N=320 

error 

error 

order 

error 

order 

error 

order 

0 

2.5292E-01 

1.9098E-01 

0.40 

1.3020E-01 

0.55 

7.9822E-02 

0.70 

LOO 

9.0170E-01 

6.8648E-01 

0.39 

4.6404E-01 

0.56 

2.8602E-01 

0.69 

1 

2.6600E-02 

4.6801E-03 

2.50 

1.0133E-03 

2.20 

2.5966E-04 

1.96 

LOO 

1.4778E-01 

3.4403E-02 

2.10 

1.1930E-02 

1.52 

3.3404E-03 

1.84 

2 

1.5883E-03 

1.8254E-04 

3.12 

2.2699E-05 

3.00 

2.8353E-06 

3.00 

LOO 

1.7729E-02 

2.7130E-03 

2.70 

3.5359E-04 

2.94 

4.4350E-05 

2.99 

3 

2.1442E-04 

1.5566E-05 

3.78 

1.0318E-06 

3.91 

6.5818E-08 

3.97 

LOO 

1.9911E-03 

2.2607E-04 

3.14 

1.5397E-05 

3.88 

9.7191E-07 

3.98 

Table  4.7 

The  KdV  equation  ut  —  3  (u'^)^  +  Uxxx  =  0.  u(x,0)  =  —2seeh^(x).  Periodie  boundary  eondition.  and  L°°  errors. 
Non-uniform  meshes  with  N  eells.  LOG  methods  with  A:  =  0,1, 2,  3.  t  =  0.5. 


k 

o 

II 

O 

00 

II 

N=160 

N=320 

error 

error 

order 

error 

order 

error 

order 

0 

2.4530E-01 

1.9004E-01 

0.37 

1.3391E-01 

0.50 

8.4650E-02 

0.66 

LOO 

1.0172E+00 

7.6826E-01 

0.40 

5.3383E-01 

0.52 

3.3672E-01 

0.66 

1 

2.7071E-02 

4.9216E-03 

2.46 

1.0581E-03 

2.21 

2.7039E-04 

1.97 

LOO 

1.4507E-01 

4.1341E-02 

1.81 

1.3916E-02 

1.57 

3.9383E-03 

1.82 

2 

2.0350E-03 

2.0344E-04 

3.32 

2.4988E-05 

3.02 

3.1228E-06 

3.00 

LOO 

2.2916E-02 

3.4702E-03 

2.72 

4.6922E-04 

2.88 

5.8972E-05 

2.99 

3 

3.2212E-04 

1.8451E-05 

4.12 

1.1715E-06 

3.97 

7.4102E-08 

3.98 

LOO 

2.8274E-03 

2.2498E-04 

3.65 

1.9437E-05 

3.53 

1.3793E-06 

3.81 

Example  4.4.  In  order  to  see  the  accuracy  of  the  method  for  nonlinear  problems  with  small  coefficient  for 
the  third  derivative  term,  we  compute  the  soliton  solution  of  the  generalized  KdV  equation  [5] 

Ut  +  Ux  +  =  0,  (4-5) 

in  — 2  <  a;  <  3,  where  we  take  e  =  0.2058  x  10“"^.  The  initial  condition  is  given  by 

u(a;,  0)  =  Asechs  (K'(a;  —  a;o))  (4.6) 

with  A  =  0.2275,  a;o  =  0.5,  and  K  =  3  j  ^  •  The  exact  solution  is 

2_ 

u{x,t)  =  ^sechs  {K{x  —  xq)  —  ujt) 

where  w  =  K  We  compute  the  solution  using  the  boundary  condition 

u{-2,t)  =  gi{t),  Ux{3,t)  =  g2{t),  Uxx{^,t)  =  gsit)  (4.7) 
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Table  4.8 

The  GKdV  equation  (4-5)  with  initial  eondition  (4-6)  and  boundary  eondition  (4-'i')-  errors.  Non-uniform 

meshes  with  N  eells.  LOG  methods  with  A:  =  0,1,2,  3.  t  =  1. 


k 

N=160 

N=320 

N=640 

N=1280 

error 

error 

order 

error 

order 

error 

order 

0 

1.6566E-02 

1.1259E-02 

0.56 

7.0817E-03 

0.67 

4.1526E-03 

0.77 

9.3056E-02 

6.6829E-02 

0.48 

4.4502E-02 

0.58 

2.7539E-02 

0.69 

1 

3.8554E-04 

6.0675E-05 

2.66 

1.1784E-05 

2.36 

2.8635E-06 

2.04 

3.2635E-03 

6.2508E-04 

2.38 

2.2689E-04 

1.47 

6.4595E-05 

1.81 

2 

8.2907E-06 

9.5348E-07 

3.12 

1.1895E-07 

3.00 

1.5290E-08 

2.96 

1.6684E-04 

2.2545E-05 

2.88 

3.0858E-06 

2.87 

3.9503E-07 

2.97 

3 

1.7005E-06 

1.3664E-07 

3.63 

3.0527E-09 

5.48 

1.9206E-10 

3.99 

1.7607E-05 

1.3291E-06 

3.72 

8.3962E-08 

3.98 

5.2861E-09 

3.99 

with  a  non-uniform  mesh.  The  result  is  contained  in  Table  4.8. 

Example  4.5.  In  this  example  we  compute  the  classical  soliton  solutions  of  the  KdV  equation 

Uf  +  ^‘^xxx  —  0-  (^•^) 

The  examples  are  those  used  in  [14]. 

The  single  soliton  case  has  the  initial  condition 

uo{x)  =  Scsech^  {k{x  —  a;o))  (4.9) 

with  c  =  0.3,  xo  =  0.5,  k  =  {\/2)^/cJe  and  e  =  5  x  10“"^.  The  solution  is  computed  in  a;  €  [0,  2]  with  periodic 
boundary  conditions,  using  elements  with  100  cells,  and  is  shown  in  Figure  4.1. 

The  double  soliton  collision  case  has  the  initial  condition 

uo{x)  =  3ci  sech^  {ki{x  —  xi))  +  3c2  sech^  (^2(2;  —  2:2))  (4-10) 

with  Cl  =  0.3,  C2  =  0.1,  2;i  =  0.4,  X2  =  0.8,  ki  =  for  i  =  1,2,  and  e  =  4.84  x  10“"^.  The  solution 

is  computed  in  2;  €  [0,2]  with  periodic  boundary  conditions,  using  elements  with  100  cells,  and  is  shown 
in  Figure  4.2. 

The  triple  soliton  splitting  case  has  the  initial  condition 

=  (4.11) 

with  e  =  10“"^.  The  solution  is  computed  in  2;  €  [0,3]  with  periodic  boundary  conditions  and  is  shown  in 
Figure  4.3. 

Example  4.6.  We  compute  in  this  example  the  KdV  zero  dispersion  limit  of  conservation  laws.  The 
equation  is  (4.8)  with  an  initial  condition 

u{x,  0)  =  2  +  0.5  sin(27r2;)  (4.12) 

for  2;  e  [0, 1]  with  periodic  boundary  conditions,  and  we  are  interested  in  the  limit  when  e  0+.  Theoretical 
and  numerical  discussions  about  this  limit  can  be  found  in  [19]  and  [23].  Here  we  are  mainly  concerned  with 
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t  =  0 


t  =  1 


Fig.  4.1.  Single  soliton  profiles.  Solutions  of  equation  (4-8)  with  initial  eondition  (4-9)  and  periodie  boundary  eonditions 
in  [0,2]  using  elements  with  100  eells.  Top  left:  solution  at  t  =  0;  top  right:  t  =  1;  bottom  left:  t  =  2;  bottom  right:  spaee 
time  graph  of  the  solution  up  to  t  =  3. 


the  capability  of  our  numerical  method  in  resolving  the  small  scale  solution  structures  in  this  limit  when  e 
is  small.  For  this  purpose  we  compute  the  solution  to  t  =  0.5  with  e  =  10“"^,  10“®,  10“®  and  10“^  using 
elements  with  300  cells  for  the  first  two  cases,  800  cells  for  the  third  case  and  1700  cells  for  the  last  case. 
We  have  verified  that  these  are  “converged”  solutions  in  the  sense  that  further  increasing  the  number  of 
cells  does  not  change  the  solutions  graphically.  These  solutions  are  shown  in  Figure  4.4.  Notice  the  physical 
“oscillations”  which  are  typical  in  such  dispersive  limits,  see,  e.g.  [19].  Clearly  our  method  is  very  suitable 
to  compute  such  solutions. 


19 


t=0  t=1 


Fig.  4.2.  Double  soliton  collision  profiles.  Solutions  of  equation  (4-8)  with  initial  condition  (4-10)  and  periodic  boundary 
conditions  in  [0,2]  using  elements  with  100  cells.  Top  left:  solution  at  t  =  0;  top  right:  t  =  1;  bottom  left:  t  =  2;  bottom 
right:  space  time  graph  of  the  solution  up  to  t  =  4. 

5.  Concluding  remarks.  We  have  designed  a  class  of  local  discontinuous  Galerkin  methods  for  solving 
KdV  type  equations  containing  third  derivatives  and  have  proven  their  stability  for  any  spatial  dimensions  for 
a  general  class  of  nonlinear  equations.  Numerical  examples  are  shown  to  illustrate  the  accuracy  and  capability 
of  the  methods,  especially  for  the  convection  dominated  cases  where  the  coefficients  of  the  third  derivative 
terms  are  small.  Efficient  implicit  time  discretizations  which  have  efficient  iterative  solvers  maintaining  the 
local  structure  of  the  method,  accuracy  enhancement  study,  and  more  numerical  experiments  with  physically 
interesting  problems  constitute  an  ongoing  project. 

Acknowledgments.  We  would  like  to  thank  Bernardo  Cockburn  for  his  valuable  help  in  the  discussion 
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t  =  0 


t  =  1 


Fig.  4.3.  Triple  soliton  splitting  profiles.  Solutions  of  equation  (4-8)  with  initial  eondition  (4-11)  and  periodie  boundary 
eonditions  in  [0,3]  using  elements  with  150  eells.  Top  left:  solution  at  t  =  0;  top  right:  t  =  1;  bottom  left:  t  =  2;  bottom 
right:  spaee  time  graph  of  the  solution  up  to  t  =  4. 


about  the  projection  S,  and  Andy  Majda  for  pointing  out  reference  [19]  and  test  cases  of  zero  dispersive 
limits  of  conservation  laws  in  Example  4.6. 
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